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1 Introduction 



Finite volume methods have been widely used in computational fluid dynamics 
for a long time; they are well adapted to the discretization of partial differen- 
tial equations under conservative form; one of their attractive features is that 
the resulting discretized equation has a clear physical interpretation [9]. In 
the framework of incompressible fluid flows, two strategies are often opposed, 
namely staggered and collocated schemes. The staggered strategy, which has 
become very popular since Patankar's book [9], remains mainly restricted to 
geometrical domains with parallel and orthogonal boundary faces. Therefore, 
for computations on complex domains with general meshes, the collocated 
strategy which consists in approximating all unknowns on the same set of 
points (called collocation points but also cell centers or simply centers), is of- 
ten preferred, even though the pressure-velocity coupling demands some cure 
for the stabilization of the well-known checkerboard pressure modes; to this 
purpose, various pressure stabilization procedures, based on improvements of 
the Momentum Interpolation Method proposed by Rhie and Chow [10], are 
frequently used [8]. 

In [3,11], a collocated finite volume scheme for incompressible flows is devel- 
oped on so called "admissible" unstructured meshes, that are meshes satisfy- 
ing the two following conditions: the straight line joining the centers of two 
adjacent control volumes is perpendicular to the common edge, and the neigh- 
boring control volumes and the associated centers are arranged in the same 
order, with respect to the common edge. Rectangular or orthogonal paral- 
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lelepipedic meshes, triangular (2D) or tetrahedral (3D) Delaunay meshes, and 
Voronoi meshes fulfill these requirements. Under this assumption, the isotropic 
diffusion fluxes can be consistently approximated by a two-point finite differ- 
ence scheme. Using this approximation for a pure diffusion problem yields a 
symmetric "M- matrix" (which ensures monotony); the stencil is limited to 
the control volume itself and its natural neighbors and it leads to the classical 
5- and 7-point schemes on rectangles and orthogonal parallelepipeds. Unfor- 
tunately, although the use of such grids considerably widens the variety of 
geometric shapes which can be gridded, it is far from solving all the critical 
needs resulting from actual problems: 

• For complex 3D domains, it is well known that the use of a large number of 
flat tetrahedra produces high discretization errors; generalized hexahedric 
meshes are often preferred: these are made of 3D elevations of quadrangular 
meshes, for which the faces of the control volumes may no longer be planar. 

• To our knowledge, there is yet no mesh software able to grid any geometrical 
shape in 3D using Voronoi or Delaunay tessellations while respecting the 
boundaries and the local refinement requirements. 

• In compressible flows, the approximation of the full tensor by the usual two 
point scheme is no longer consistent even on admissible meshes; multi-point 
approximations are therefore required. 

• Boundary layers are classically meshed with refined grids, so that the dis- 
cretization scheme should be able to deal with non-conforming meshes. 

Whereas there is no real difficulty to discretize the convective terms for general 
non-conforming grids, writing accurate diffusion approximations, particularly 
relevant for low Reynolds (Peclet) flows, is still a challenge on such meshes. 
In the early 80's, Kershaw [7] first proposed a nine-point scheme on structured 
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quadrilateral grids by using the restrictive assumption of a smooth mapping 
between the logical mesh and the spatial coordinates. Since then, numerous 
works have been published to efficiently solve the diffusion equations in general 
geometry (see [1] for a review of recent papers). The drawbacks of the actual 
schemes for diffusion are often linked with one or several of these key points: 

• a non-local stencil (quite dense matrices); 

• cell-centered but also face-centered unknowns (large matrices); 

• non-symmetric definite positive matrices (loss of the energy balance); 

• loss of the convergence or of the accuracy on some particular grids; 

• loss of monotony for solutions in purely diffusive problems (the resulting 
matrix is not an "M-matrix" ) . 

We focus in this paper on the approximation of the Navier-Stokes and energy 
equations under the Boussinesq assumption, using a new scheme for diffusion 
terms. This scheme is shown to provide a cell-centered approximation with 
a quite reduced stencil, leading to symmetric definite positive matrices and 
allowing a mathematical proof of convergence. Although the diffusion matrix 
may not be shown to be an M-matrix in the general case, the maximum prin- 
ciple is nevertheless preserved in our numerical three-dimensional simulations. 
In this scheme, the discrete pressure gradient and the non-linear contributions 
are approximated so that the discrete kinetic and energy balances mimic their 
continuous counterparts. Indeed, the pressure gradient is chosen as the dual 
operator of the discrete divergence, and the discretization is such that there is 
no contribution of the non-linear velocity transport in the increase of kinetic 
energy. In order to suppress the pressure checkerboard modes, the mass bal- 
ance is stabilized by a pressure term which only redistributes the fluid mass 
within subsets of control volumes, the characteristic size of which is two or 
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three times the local mesh size. 

The remainder of this paper is divided into four sections. In Section 2, the con- 
tinuous formulation is presented in the framework of free convection. Section 3 
presents the discrete scheme for general non-conforming meshes, with illustra- 
tions in the simplified case of uniform rectangular grids, and some mathemat- 
ical properties. The fourth section is devoted to the numerical validation, first 
with analytical solutions and then with a classical natural convection problem. 

2 Continuous formulation 

Let d be the dimension of the space (d — 2 or 3) and let Q C M. d be an 
open polygonal connected domain. For x e Q, our aim is to compute an 
approximation of the velocity u(x) = Y^i=i u ( x ) e ij the pressure p(x) and 
the temperature T(x), solution of the steady and dimensionless Navier-Stokes 
and energy equations under the Boussinesq approximation: 

-PrAu + Vp + (u ■ V)« - RaPrTe 3 = f(x) in Q, (la) 

—AT + (u ■ V)T = g(x) in fi, (lb) 

divu = in f2, (lc) 

where e-j indicates the vertical upward direction, f(x) = Y,f=i f { x ) e i an d 
g(x) are dimensionless regular functions modeling source or sink in the mo- 
mentum or heat balances; Pr and Ra denote the Prandtl and Rayleigh num- 
bers respectively. We consider the case of the homogeneous Dirichlet boundary 
conditions for the velocity and of the mixed Dirichlet-Neumann boundary con- 
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ditions for the temperature. These boundary conditions read as follows: 
u(x) = a? G T, 



< T(x) = T b (x) 



x e r l5 



(2) 



-VT(as) • n(x) = q b (x) x G T 2 



where Ti,]^ are subsets of the boundary T of the domain Q such that T\ fl 
T 2 = and r\ U T 2 = T, and n(x) is the outward unit normal vector to the 
boundary. We assume that T b is the trace on Ti of a function again denoted 
T b such that Tb G H l (Q), and we define the functional space (f2) = {T G 
H 1 (Q);T(x) = on Ti}. Then a weak formulation of equations (la-lc) with 
boundary conditions (2) reads: 

Find u G HliVt) d , p G L 2 (Q) with J n p(a;)da; = 0, and T with T - T b G 
ifp x (f2) such that 



Pr / Vu : Vvdx — / pdivvda3+ / div(u <g> u) ■ vdx 
Jn Jn Jn 



(3a) 



-RaPr / Te 3 • vdx = [ fix) ■ vdx, \/v G H^(fl) d 
Jn Jn 



[ VT • VOdx + / dw(uT)6dx 
Jn Jn 



g(x)9dx- [ q b (x)6(x)dx, V0 G Hi (£l). 
n Jr 2 11 

divu(x) = for a.e. x G fl, 



(3b) 



(3c) 



Although our discretization scheme belongs to the finite volume family, we 
shall also be using the weak form (3a-3c) in our discretization. Indeed, the dis- 
cretization of the diffusive terms —PrAu in (la) and — AT in (lb) is obtained 
by the construction of a discrete gradient which is then used to approximate 
the term Pr J u Vu : Vvdx in (3a) and J n VT ■ V6dx in (3b). 
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3 Numerical scheme 

In this section we present the discretization scheme for Problem (l)-(2) under 
its weak form (3). The next paragraph is devoted to the notations for general 
discretization meshes and to the description of the discrete degrees of free- 
dom. We then describe the approximation of the diffusive terms (Sec. 3.2). 
Because of the collocated choice of the unknowns, a stabilization is needed. 
The stabilization we choose is imposed on the mass flux (rather than the over- 
all balance) and also appears in the momentum and energy equations through 
the convective contributions: this is described in Section 3.3. It also involves 
the choice of some coefficients which are defined in Section 3.4. The complete 
discrete problem is finally given in Section 3.5, and some of its mathematical 
properties sketched in Section 3.6. 

3.1 Mesh and discrete spaces 

We denote by T> = (Ai,£,V) a space discretization, where (see Fig. 1): 

• M. is a finite family of "control volumes", i.e. non empty connected open 
disjoint subsets of Q such that fl = Uk^mK- For any control volume K £ 
A4, we denote by dK = K \ K its boundary, mx > its measure (area if 
d = 2, volume if d = 3) and hx its diameter (that is the largest distance 
between any two points of K). 

• £ is a finite family "edges" (d = 2) or "faces" (d = 3) of the mesh; these are 
assumed to be non empty open disjoint subsets of Q, which are included in 
a straight line if d = 2 or in a plane if d — 3, and with non zero measure. 
We assume that, for all K £ Ai, there exists a subset 8k of £ such that 
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OK = \J a ££ K ~a. The set £ is assumed to be partitioned into external and 
interior edges (d = 2) or faces (d — 3): £ = £ mt U £ cxt , with a C dQ for any 
a G £ cxt and cr C f2 \ <9f2 for any a G £i nt . Any boundary edge a is assumed 
to belong to a set £k for one and only one K E A4; any interior edge a 
is assumed to belong to exactly two sets £k and El with K ^ L, and in 
this case a is included in the common boundary of K and L, denoted K/L. 
Note that there are cases in which Kj L includes two or more edges or faces, 
see for instance the third mesh for the unit cube, section 4.1, and Figure 
5b. We also assume that, if a G £ cxt , then either a C Y x or a C T 2 . For all 
er G £, we denote by and m^ the barycenter and the measure of a. For 
all K G A4 and a G £k, we denote by the unit vector normal to a 
outward to K. 

• V is a family of collocation points V = (x K ) KeM of Q which is chosen such 
that for all K G M. and for all x G K, the property [xk, x] G K holds. 
Note that this choice is possible for quite general polygons, including those 
with re-entrant corners, see Fig. 1. The Euclidean distance dx, a between Xk 
and the hyperplane including a is thus positive. We also denote by Ck,u the 
cone with vertex Xk and basis a. 

Next, for any a G £; nt , we choose some real coefficients such that the 

barycenter x a of a is expressed by 

** = E E /?. L = i. (4) 

In three space dimensions, it is always possible to restrict the number of 
nonzero coefficients (3% to four (in practice, the scheme has been shown to 
be robust with respect to the choice of these four control volumes, taken close 
enough to the considered edge). 
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Note that in the case of an uniform rectangular grid which is depicted in 
Figure 2, an obvious choice for the coefficients (3% is obtained by noticing 
that x i+ i/ 2 j = (xjj+Xj+ij)/^ and x i>j+1/2 = (x i)i + x i j + i)/2. Thus for any 
edge a, only two coefficients need to be nonzero. 

We now define the finite dimensional space M. M x M £ (where an element v G 
WL M x ~R £ is defined by the family of real values ((vk)keMi ( v a)aee)) and the 
following subspaces: 

• X v = {u G R M x R £ , Vcr G S int ,u a = El&m PSl} (the dimension of X v 
is the number of control volumes plus that of boundary edges), 

• X"q = |m G X 25 , Vcr G £ cxt , ^cr = o| (the dimension of X® is the number of 
control volumes), 

• x£ lfi = [9 G Vcr g £ cxt n ri,0 ff = 0} (the dimension of Xj? is the 
number of control volumes plus that of boundary edges on 1^). 

3.2 Discretization of diffusive terms 

Let us first define a discrete gradient for the elements of X v on cell K G M. 
We set, for any u G X v and K G M: 

V K u = — m <r( u o - u K )rt K)a . (5) 

Note that this is a centered gradient. 

As an illustration, consider the case of the two dimensional uniform rect- 
angular grid depicted in Fig. 2, let u G X v , and choose the natural choice 
(3% = 1/2 if a is a side of L and otherwise. Then with the (natural) 
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notations of Fig. 2, one has: 



1 1 ^ 

— (Ui+ij -Ui_ij) 

z/ we app/?/ t/izs formula to the element 1^ ^ G u>z'i/i a// o/ components 
equal to zero except for the one associated with K it j which is equal to 1, we 
get that the vector V^Ir-^ is zero for all control volumes L except those 
neighboring Kij, as shown on Fig. 3a. Considering the checkerboard so- 
lutions on uniform rectangular grids, the above expression shows that this 
discrete gradient may vanish for some non- constant functions. 

An approximation of the diffusive terms Fr f n Vu : Vvda: in (3a) and f n VT- 
V^dx in (3b) using this discrete gradient (5) would yield a non-coercive form. 
Thus, we shall work with a modified gradient, defined in (6)-(7) below. 

To this end, for all a G £k we first define Rk,o u £ K which may be seen as a 
consistency error on the normal flux, by: 

\fd 

Rk,cjU = -j \U a ~U K — VjfM ■ (x a - x K )) . 

One may note that Rk^u = if uk and u a are the exact values of a linear 
function at points x K and x a , for all K and a G £k- We then give the following 
expression for a stabilized discrete gradient of u G X v in each cone CV jCT : 

Vjf, ff M = V K u + Rk,c,u n K ,a- (6) 

In the case of the uniform rectangular grid given in Fig. 2 and for K = K it j 
and a = erj+i^j, we find: 



10 



The stabilized discrete gradient (6) applied to the above defined element 1k u 
of X v provides nonzero contributions on the triangular subcells Cx,a of the 
cell Ki j and of its neighbors (Fig. 3b). 

The global discrete gradient is chosen as the function Vpti: 

V v u(x) = V K ,crU, for a.e. x G C K>CT , WK G M, Va G £ K . (7) 

We then plan to approximate the term f n Vu(cc) • 'Vv(x)dx by: 

/ V v u(x)-V v v(x)dx = £ E ^^V K>( ,u-V K>a v, Vu,veX v . (8) 

In fact, it is shown in [4,5] that this expression defines a symmetric inner prod- 
uct on X v , and provides a good approximation for / n Vu(a;) ■ Vv(x)dx; this 
approximation may be seen as a low degree discontinuous Galerkin method. If 
one seeks a finite volume interpretation of this scheme, it is possible, expressing 
Ufj and v a for all a G £; n t thanks to the relations (4), to show that 

/ 'Vt)u(x) ■ Vvv(x)dx = 

x (9) 

where for any K G M, Mk is the subset of cells playing a part in the barycenter 
expression of x a , for all edges of the cell K and of the neighbors of K, i.e. 
Mk = {M G M\ffi ^ 0,W G £ L ,VL e M a ,W G £ K }; F K , L (u) is a linear 
function of the unknowns (ul)l£M which is such that Fk,l(u) = —Fl^k{u). 
In the general case, the expression of Fk,l(u) is rather complicated (see [5]). 

In the case of the uniform rectangular grid of Fig. 3, this expression simpli- 
fies into the usual two point flux; for instance, the flux from K^j to -f^+ij 
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reads: 

F KiJ ,K i+1 Au) = h 



"y ^ • 



More generally, a two point flux is also obtained in the two or three di- 
mensional non uniform rectangular cases. Indeed, locating x^ at the cen- 
ter of gravity of the cell K, the relation {x a — XK)jdK,a = tik,(t holds. 
It is then possible to write x ff = (di,a^K + dx^i) / {dL,a + dx,a) and 
u a {v) = (dL tG u K + dK^ u L)l{dL,a + dx,a) f or all °~ such that a C K/L, 
and for all u G X v . Using the identity 

m CT (£c CT - x K )n^ i(T = m K l 

where 1 designates the transposition and I the identity matrix, we obtain [5, 
Lemma 2.1]: 

Vt>u(x) ■ Vt>v (x)dx = Yl 1 ri — i u L ~ u K )(v L - V K ) 

A ae& mU *CK/L d K,« + d L,a 

+ E -j-^-iUa - U K )(v a - V K ). 

Then the previous relation leads to define Mk as simply the set of the natural 
neighbors of K , and to define the fluxes by the natural two-point difference 
scheme, in the same manner as in [3,11]: 



xxi 

Fk,l( u ) = 1 — > — ( u k ~ u l) for a e £ int , o C K/L 

XXX 

F K ,*(u) = 1 — ( u k Ucr) for a e S cxt PI S K . 



(10) 



The classical and cheap 5- and 1-point schemes on rectangular or orthogonal 
parallelepipedic meshes is then recovered. An advantage can then be taken 
from this property, by using meshes which consist in orthogonal parallelepi- 
pedic control volumes in the main part of the interior of the domain, as 
illustrated by the cone-shaped cavity (Fig. 5c) in Section 4-1- 
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Note that the approximation of — J K Audx is obtained by letting vk — 1, 
Vl = for L 7^ /<" and u CT = for a G £ cx t in (9): 

- f Audx ~ ^ F A - L ( M ) + ^ F Kja (u), 
K LeN K cre£ K n£ C xt 

so that we may define an approximate Laplace operator A-p by the constant 
values A^u on the cells K: 

-A K u = — f E f k,l{u) + E ^av(m) ] . (11) 

mK \LeAf K <re£ K n£ c: 



- cxt 



In t/ie case of the uniform rectangular grid of Fig. 3, this discrete Laplacian 
leads to the usual five point formula: 

-A Ki j U = —2 (2Uij - U i+ xj - Ui-xj) + (2u id - U id+1 - MjJ-l) . 

ll X fly 

In the general case, the stencil of the discrete operator on cell K is defined by 
Mk (see relation (11)) and therefore depends on the way the barycenters x a 
are computed. For general grids, the equation for a given cell usually concerns 
the unknowns associated to itself, its neighbors, the neighbors of its neighbors 
and possibly some additional adjacent cells. The resulting matrix is usually not 
an "M- matrix" , except on particular meshes such as conforming orthogonal 
parallelepipeds, in which case we obtain the usual two point flux scheme, as 
previously pointed out. 



3.3 Pressure-velocity coupling, mass balance and convective contributions 



For all u G (X^) d , we define a discrete divergence operator by: 



div K u = — m a u a ■ riKa, Vif G M. (12) 
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where u a = J2i=i u a^ e i- Notice that 

div x u = £(V^)« 

i=l 

with Vif«w defined by (5). 

In the case of the uniform rectangular mesh of Fig. 2, this operator reads: 



di\ Kij u = 




We then define the function div^u by the relation 

divx>u(x) = divert, for a.e. x G K, WK G M.. 

The discrete gradient operator used for the pressure gradient is defined as the 
dual of this divergence operator. More precisely, we mimic at the discrete level 
the (formal) equality f n p divv dx = — f n Vp ■ v dx. The discrete equivalent 
of f n p divv dx reads J2LeM m L Pl div^i; with v G (X v ) d and p G X v ; we 
then define the discrete pressure gradient which we denote V ' kP (on cell K) , 
such that 

J2 rn L V L p ■ v L = - J2 ™- L p L di\ L v, V v G (X v ) d . 
LeM LeM 

^From the definition of the divergence (12) and of X v , we thus seek (V lP)l^m 
such that: 

^2 m L V L p-V L = - ^ p L m ^ Yl Pa V M- n L,a- (13) 

LeM LeM a-eS L n£ iBt MeM 
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Taking for v the element of (X v ) d with components v[ — 1 if j — i and 
L = K, and otherwise , we thus get: 

m K V K p = - J2 PL J2 m <Tpa n L,a 

= Yl m a (3^(p M - PL)n L>a . (14) 

a€S inU aCL/M 

Remark that 'VkP is neither constructed with the discrete gradient (5) nor 
with the stabilized one (6); its expression is the dual form of the divergence 
(12). 



In the case of non uniform rectangles (d = 2) or parallelepipeds (d = 3) 
with collocated points at the gravity center of the cells, and for all a G £ in t 
with a C K/L, we only need two non-zero coefficients (3^ : (3^ - 



dn,cr+d.L,, 

— TYt vrefrtve rel nt') nrt i 1 /. I TpAnrpQ in 

d.K,ir+dL,, 



and (3% = , dK .^ — . Therefore, relation (14) reduces to 



va K V K p= ~T1 — (PL - PK)n K>(T . 



Using p K J2aes K m CT n^ i(T = 0, 



^ x - d L ^p L + d K a p K x 
m K V K p = 2^ HV — — n K ,a + 2^ ra a p K n K , a 

ae£ K ,aCK/L + dL ,° a££ K n£ e: 



-'CXt 



For a uniform grid and a control volume without boundary faces, the above 
expression resumes to 

e \- pl + pk 

m K V K p = 2^ m CT % iffI 

a££ K ,aCK/L 

which provides, in the particular case of Figure 2, the usual formulation 



h x h y V K p = hJ''- Ll J" + hj''' 1 ' l ''' J ' 
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As recalled in the introduction of this paper, a pressure stabilization method 
is implemented in the mass conservation equation in order to prevent from 
oscillations of the pressure, as for instance in [2] in the finite element setting, 
[8,10] in the finite volume setting. The originality of our approach is that we 
directly include the stabilizing diffusive pressure flux in the approximated mass 
flux, so that it will appear not only (as usual) in the mass equation, but also 
in the momentum equation through the non-linear convective term. From the 
mathematical point of view, this helps in obtaining simple estimates on the 
velocity and pressure, but more importantly, it ensures that the contribution 
of the discrete non-linear convective term to the kinetic (and thermal) energy 
balance is zero, just as in the continuous case. Let us define the stabilized 
mass flux across a C K/L by 



where (A CT )o- e £ int is a given family of positive real numbers, the choice of which 
is discussed below. Note that the quantity X a (px — Pl) ma y be seen as a 
numerical pressure diffusion flux, and that the overall numerical flux remains 
conservative, that is, if a C K/L then $^ CT (u,p) + ^ a (u,p) = 0. 

In the case of the uniform rectangular mesh of Fig. 2, the expression of this 
flux through a vertical edge o" i+1 / 2 j reads: 



We then use the modified flux, in order to define a stabilized centered transport 
operator which is defined, for all u e (X^) d , w G X v and K e Ai, by 




(15) 





1(3 



An interesting remark is that, in the case where the mass balance equation in 
the control volume K is satisfied, that is: 

div^(l,«,p) = -J- E ^ <a (u,p) = 0, 

then 

E ®K, a (u,p)w K = 0, 
so that the following relation also holds: 

ae£ K ,<TCK/L Z 

We shall use this latter form in the practical implementation, in particular 
in the discretization of the non-linear convection term. Indeed, it is more 
efficient when computing the Jacobian matrix of the momentum equation, 
since it avoids summing up values of the same amplitude. 

In the particular case of the uniform rectangular grid of Fig. 2, the sum- 
mation Yjae£ K ,acK/L involves the four edges between the control volume 
K = Kij and its neighbors L = Ki+ij, Ki-ij, Kij + i, Kij-i. If L = i^i-ij 
the expression which appears in the summation reads: 

h y I 2 M-i/2,j[Pi,j -Pi-i,j) I 2 • 

When the local grid Reynolds (or Peclet) number is much larger than 1, an 
upwind scheme must be applied that consists in substituting div^(w, u,p) by 

div^ up O, u,p) = 

— (max(^ <a (u,p),0)w K + mm(<l> x Ka (u,p),0)w L ). 

mR a££ K ,aCK/L 

In both the centered or upwind cases, the functions div^(w, u,p) and div^,' up (w, 
are defined by their constant values in each control volume. For u,w e (X^) d , 
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we also define the centered vector divergence operator div v (w, u,p) such 
that the i-th component of dW^(w, u,p) is equal to div^,(u;M, u,p), for i = 
1 , . . . , d; a similar expression applies for the upwind vector divergence operator 
div^,' up '(w,u,p). 

3.4 Choice for the parameters (A CT ) CT g£: int 

Different strategies can be applied to define the parameters (X a ) ae £ int . Amongst 
all of them we applied the "cluster stabilization method" [3,11]; it consists in 
constructing a partition of A4, denoted Q, and setting A CT = A > if there 
exists G G Q (such G C M. is called a cluster) with a C K/L, K and L 
belonging to G, and A CT = otherwise. Here is an example of an algorithm 
creating a cluster partition: 

(1) for all cells K e Ai, initialize a new cluster if K and its neighboring cells 
do not already belong to a cluster; 

(2) for any remaining isolated cell L, connect it to the closest cluster having 
the largest number of common edges with L. 

This algorithm is now applied to the mesh of figure (4a), the cells being de- 
scribed from left to right, from the lower to the upper row. Figures (4b) and 
(4c) illustrate the first cluster and the set of the clusters at the end of the first 
step of the algorithm. Figure 4d shows the clusters after the isolated (unnum- 
bered) cells of figure 4c have been connected. The choice of the value A depends 
on the physical problem and it must be chosen both large to avoid the ap- 
pearance of spurious modes and small to preserve an accurate approximation 
of the mass equation. In [11], comparisons with other stabilization methods 
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were performed and the results showed that solutions are little sensitive to 
the value of A. For the natural convection example presented in section 4.2, 
A = 1(T 8 . 



3.5 Resulting discrete equations 



We denote by T^x> the element T £ X v such that Tk = for all K £ Ai, 
T a = for all a £ E int and all a £ £ cxt with o C T 2 , and, for all a £ £ ext with 

T a = —l T b (x)ds(x). ( 16 ) 



m CT Jo- 
Let ?%vt(f2) C L 2 (Q) denote the set of functions which are constant in each 
K £ Ai] for any function q £ T^^(fi), we shall denote by qx its constant 
value on K £ Ai. We then define the mapping Pm : X 15 — > H.m(Q) by 
w £ X v i — ► Pv(t) with P M u(a;) = u_k for a.e. x £ K and all if £ .M. We also 
define the mapping Ps : — > L 2 (T) by u £ X v i— > Pgv with Psv(x) = v a 
for a.e. x £ ex and all a £ £ ext . 

Let us then use the previously defined discrete operators to formulate a discrete 
approximation to problem (3): 

Find u = (u») i=M £ (X^) d , p £ H M (n) with J n p(x)dx = Ek^m ™kPk = 
and T — £ Xj? j0 such that: 



Pr / VpW : Vp?} dx — p divp-u dx + I div^,(it, it,p) ■ Pm v dx 

Jn Jn Jn ^\ 

-RaPr / P M T e 3 ■ P M v dx = [ f ■ P M v dx, \/v £ (Xg) d , 
Jn Jn 
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/ V V T ■ V v dx+ f div£(T, u,p)P M 6 dx 
Jn Jn 

= ( g p M e dx - [ qb p £ e ds, \/e e x£ , 

J a Jt 2 



div£(l,u,p) = a.e. in Q. ( 19 ) 

We then deduce from (17) the d discrete momentum balances over the control 
volume K, letting v^' = 1 in K, and otherwise; these equations read, in 
vector form: 



-Fr niA'AjfM + Y ™-aPa (j>M ~Pl) n L ,a 

ae£ int ,crCM/L 



(20) 

+ Y *KA u >P) UK t UL " ***** m KT K e 3 = [fdx 

(where —A^u is the vector valued discrete Laplace operator defined by (11) 
for each of its components). Similarly, we deduce from (18) the discrete energy 
balance over the control volume K, letting 9 = 1 in K, and otherwise; this 
equation reads: 



-m K A K T+ Y ^,Au,p) J ^± = j K gdx. (21) 



T K + T L 

a££ K ,aCK/L 

Recall that, for all K G Ad, and all a G such that a C Y\, the Dirichlet 
boundary condition (16) is given. We deduce from (18) the relation imposed 
by the Neumann boundary condition for the thermal flux, letting 9 a = 1 and 
otherwise, for some a G £k with cr C ^ 



FkAT) = I q b {x)ds{x). 

J a 



[22] 
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Note that the above relation is natural, accounting for the fact that Fk,o{T) 
approximates the heat flux at the edge a. Finally, we write (19) in a given 
control volume K: 

E **>'P) = °- (23) 

ae£ K ri£ int 

The stencil of the scheme (20-23) is always determined by that of the diffusion 
operator. 

As previously mentioned, in the case of orthogonal quadrilateral or par- 
allelepiped grids, the diffusion flux through an edge o gives the classical 
two-point difference scheme (10) and the bary center coordinates are simply 
calculated by an arithmetic average, x a = {d Lcj x K + d K,a x L)l (d>L,a + dx,a)- 
Hence, in this case, equations (20) -(23) then read, for a given K e Ai: 

Eat + E i K VT = J R gd 
M (K) - 



with: 

cr€£ K ,aCK/L Ct ^' <T + CtL > <7 

n (K) v / d L>a u K + djc^UL w A UK + U L 

*e£ K ,aCK/L \ d K,a + d Lia ) 2 

= — RaPr m K T K e 3 , 



x. 



E A t= E T^T" + E m -— 3 

ae£ K ,aCK/L ~T L,cr ae£ K n£ eKt UK ' a 



UK) _ V- Tr , f ^.^AT + 4, g «L » , A T K + T L 

E U -VT - 2^ m <M j T~l n *> + - PL) 7, : 

*e£ K ,*CK/L V d K,a + d L<a ) 2 
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m <r n Ki<T 



d L ,gU K + d K ^u L 



+ K(j>K ~Pl) 



) 



and with the Dirichlet boundary (16) applied to a G Ti and the Neu- 
mann boundary condition on a C T 2 being reduced to m a (T a — T K )/dx,a = 

~ IaQb(x)ds(x). 

3. 6 Some mathematical properties 

The system of discrete equations (20-23) is a system of non-linear equations. 
The mathematical proof of the existence of at least one solution can be shown 
in the particular case T& = and g& = 0, which we consider in this section. 
Indeed, in this case, we can show some a priori bounds on T and u. We first 
let 9 = T in (18). Using the relation 



Thanks to a discrete Poincare inequality which follows from [5, Lemma 5.3], 
we get that there exists Ct, only depending on the regularity of the mesh and 
on g, but not on the size of the mesh, such that 




which results from (19), we get 




We then let v = u in (17). We get, thanks to (15) and (19) 



Pr||V©u||( L 2 (n)d) d + 



m a X a (p L -p K ) 2 



ae£ int ,<TCK/L 




22 



Again using the Poincare inequality, we conclude that there exists C u , only 
depending on the regularity of the mesh, on Ra, Pr, / and g, but not on the 
size of the mesh, such that 

W^vu\\( L 2( n -)dy < C u . 

Hence, using the topological degree method, we can prove the existence of 
at least one solution. Moreover, these inequalities are then sufficient to get 
compactness properties, which show that, from a sequence of discrete solutions 
with the space step tending to zero, we can extract a converging subsequence, 
for suitable norms. Then we can prove that the limit of this subsequence has a 
sufficient regularity, in relation with the weak sense provided by (3). It is then 
possible to pass to the limit on (18), (17) and (19), using test functions which 
are interpolation of regular ones. We then get that the limit of the converging 
subsequence satisfies (3). 

4 Numerical validation 

Numerical implementation is performed for three dimensional domains. We 
first validate our results on known analytical solutions which allow us to com- 
pute the scheme's order of convergence. We then turn to a natural convection 
case which is referenced in the literature. 

The domains considered are either the unit cube or a circular centered cone. 
In both cases, we use structured (rectangular and non rectangular) meshes, 
and we denote by N the number of cells in each of the three space directions. 
The set of non-linear equations (20-23) is solved by an under-relaxed Newton 
method where the unknowns are the velocity %, the pressure px arid the 



23 



temperature Tk and T a for all K e A4 and a e £ ex t H The solutions of the 
linear systems are computed with a parallel Generalized Minimal RESidual 
method provided by the scalable linear solvers package HYPRE with a pre- 
conditioning based on the block Jacobi iLU factorization carried out by the 
Euclid library [6] . 



4-1 Analytical solutions 



We consider two closed cavities, cubic or cone-shaped, in which the fluid flow 
and the heat transfer are known a priori. Let p re f , u re f and T re f be some known 
pressure, divergence free velocity and temperature fields; we then compute / 
and g by Eqs. (la, lb) where we have set u = u re f, p = p rc f and T = T re f . 
For any regular function ip (ip = [u^\ i = 1, • • • , 3, p or T), the scheme's 
relative accuracy for the usual L°°, L? and H l norms is measured by 

maxxex H>K ~ ipref(x K )\ 



e 2 W 



1 /2 

' T,K&M T0a- K \V K xjj - V^rcf(^^)| 2 \ 

E^ e ^( m ^lWref(a;K)| 2 / 



(24) 



where | • | denotes the usual Euclidean inner product in M 3 . For each of the 
above relative error, the scheme's order of convergence is defined by the mean 
slope of the logarithm of the relative error as a function of the logarithm of 
the largest cell diameter max^- gj v( hx, the slope being calculated by a least 
square method. 

Three different meshes are studied for the unit cubic enclosure. Except for the 
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last mesh, the points Xk are located at the gravity center of the cells. 

(1) The first mesh is an uniform mesh consisting of orthogonal parallelepipeds 
of size 1/iV 3 . 

(2) The second one (Fig. 5a) is constructed by a smooth mapping between the 
uniform mesh and the spatial coordinates [1]. The vertices x s (i,j,k) = 
ixjp(i,j,k)) of the elementary distorted cubes are defined by: 

V J 1=1,- ,3 

V{i,j,k) e N([l,m + 1]) x N([l,n 2 + 1]) x N([l,n 3 + 1]), 
xW( i>i> A ; ) = i- o B (^^) 

-P<U*) = ^ + 0.1 sin (^) sm (^il) 
^(,^) = ^ + 0.1sin(^) S in(^). 

(3) The third mesh (Fig. 5b) is constructed from the first one in the follow- 
ing way: the points Xk remains at the gravity centers of the basic mesh 
whereas the vertices x s of the cells are randomly displaced in each space 
direction at most of 0.45/ N. Unlike the previous mesh, which consists in 
hexahedra with plane faces, the four edges of a face are now not included 
in a same plane, with the exception of edges which belong to the bound- 
aries of the cubic domain. Since the consistency of the discrete gradient 
defined in (5) (and therefore of that defined in (7)) holds if the faces a 
are plane, we replace each of the non planar faces by two triangular faces. 

The second enclosure, the cone-shaped cavity, is bounded by the lateral surface 
((xW - 0.5) 2 + {xW - 0.5) 2 = ((6 - 5a^)/12) 2 for x^ G [0, 1] and by two 
plane discs (x (1) - 0.5) 2 + (x (2) - 0.5) 2 < 1/4 for x (3) = and (x (1) - 0.5) 2 + 
(x^ — 0.5) 2 < 1/12 2 for x^ = 1. To make the most of the scheme which 
locally reduces into a 7-point difference scheme on parallelepiped cells, the 
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mesh consists of cubes which were cut to match the lateral curved boundary. 
Thus, the mesh error tends quadratically to zero with respect to the mesh 
size. To avoid too large differences of volume sizes between adjacent cells that 
may deteriorate the numerical accuracy, the boundary cells having a volume 
less than 0.1 /{rii) d are merged into adjacent cells. 

We are first interested in the Poisson problem for a scalar variable (Eq. (21) 
with u = 0) where r\ = V and T 2 = 0. It was first checked that the er- 
rors obtained with a linear analytical solution on the different meshes and 
cavities are of the order of the computer accuracy, even for the coarsest 
grids. The next analytical test consists in choosing the reference solution 
T rc f (x^\ x^ 2 \ x^) = sin(7rx^ 1 ' ) ) cos(nx^) cos(7rx®) with appropriate Dirich- 
let boundary conditions (Fig. 6a-c). The orders of convergence are reported 
in table (1). The accuracy of the scheme is close to 2 when considering the 
L 2 -norm and it slightly decreases with the L°°-norm but always remains larger 
than 1.50. The order of convergence for the gradients (if 1 -norm) is larger than 
I. 

Next, we examine the convergence behavior of the isothermal Navier-Stokes 
equations by setting u ref (:r) = V A £f =1 (4x (1) (x (1) - l)) 3 (4x (2) (x (2) - l)) 4 
(4x^ 3 \x^ — l)) 5 ej and Pref(^) = cos^x^) cos(7ra;( 2 )) cos("7ra( 3 )) in (Eq. la) 
with Pr = 1 and Ra = (note that the dimensionless writing of the equa- 
tions is meaningless because the current reference velocity is related to the 
thermal diffusivity which has no reason to appear for isothermal problems. 
Another velocity reference should be used, based on the viscous diffusivity). 
Table (2) indicates that the convergence rates of the velocity components are 
larger than 1.90 on the three finer meshes when the relative error is based on 
the L 2 -norm and first order accurate for the pressure for distorted meshes. In 
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accordance with the diffusion problem when the L°°-norm is used, the orders 
of convergence slightly decrease for the velocity but a convergence rate larger 
than 1.60 is still measured. The convergence rates of the gradients are better 
than the expected first order. Unsurprisingly, the L°° and if 1 -norms of the 
pressure do not tend to zero with the mesh size because it simply appears in 
the momentum equation as Lagrangian multiplier of the mass equation. Thus 
the only guaranteed convergence for the pressure is based on the L 2 -norm. 

4-2 Natural convection problem 

We consider an air filled unit-cubic enclosure with isolated walls except the two 
face to face vertical isothermal surfaces at = and 1. The governing fluid 
flow equations are solution of system (1) and (2) with fix) = 0, g(x) = 0, 
q b (x) = 0, T(0,x (2) ,x (3) ) = -0.5 and T(l, x (2) , x (3) ) = 0.5. The Prandtl and 
Rayleigh numbers are fixed to Fr = 0.71 and Ra = 10 7 and the stabilization 
parameter is chosen equal to A = 10~ 8 in the mass equation. Because very 
small boundary layers take place along the walls, the vertices are located at 
the Gauss-Lobatto points and the collocation points xk at the gravity centers 
of the cells. To also study the effect of non-cubic meshes, the coordinates of 
the previous defined vertices are randomly moved, for each direction, of a 
magnitude at most equal to 0.45 the size of the cell in this direction. 

Table (5) presents the maxima of the velocity components, the average Nusselt 
number on the isothermal walls and their relative differences with respect to 
reference data [12]. For the cubic meshes, the relative differences with respect 
to reference data seem to tend to zero: for the finer mesh, our results depart 
from less than 1% with respect to the reference values. Although the solutions 
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are less accurate for the shaken meshes, a convergent behavior is also observed. 
The relatively large differences measured on u^}, in comparison with the other 
velocity components, are probably the result of the flow shape in which the 
main motion occurs in the (ei, ea)-planes with a small secondary flow in the 
transverse planes. It is also interesting to note the good accuracy of the average 
Nusselt numbers. This accuracy is essentially obtained thanks to the use of the 
stabilized mass flux ^^ a (u,p) in the heat transport expression which ensures 
the conservation of the average heat flux balance on the boundaries of the 
cavity. 

5 Conclusion 

In this paper we presented a new scheme which is well suited for the simu- 
lation of incompressible viscous flows on irregular and non-conforming grids. 
This possibility seems to open a large field of new applications (grid refinement 
as a function of an a posteriori error computation, free boundaries, . . . ). We 
emphasize that the convergence of the scheme may be proven mathematically, 
and that the obtained numerical results are accurate. Although we presented 
this scheme in the steady case, its extension to transient regimes is straight- 
forward. In this latter case, one should consider optimizing the linear solving 
step by using suitable projection algorithms. 
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